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1.0  INTRODUCTION 

The  overall  goal  of  this  project  is  to  develop  an  experimentally  validated 
computational  model  of  the  eye  and  apply  the  model  to  evaluate  the  stresses  and 
deformations  incurred  by  the  eye-wall  and  critical  ocular  components  from  blast 
overpressures,  and  to  investigate  the  interaction  between  the  standard  issue  eye 
armor  and  the  blast  wave,  and  its  effect  on  the  mechanical  loading  of  the  eye.  The 
model  will  be  developed  based  on  the  following  working  hypotheses.  1)  The 
anisotropic  mechanical  properties  of  the  cornea  and  sclera  as  derived  from  the 
collagen  structure,  are  critical  to  modeling  the  interaction  of  the  blast  wave  and  the 
globe.  2)  The  mechanical  behavior  of  the  cornea  and  sclera  under  dynamic  (high- 
rate)  loading  is  significantly  different  than  the  under  qua-sistatic  (slow-rate) 
loading.  3)  The  surrounding  environment  of  the  globe,  including  the  extraocular 
tissues  of  the  orbit  and  bony  facial  features  is  important  in  determining  the  effects 
of  blast  loading  on  the  eye.  In  this  second  year  of  the  project,  we  focused  on 
developing  a  computational  model  to  assess  the  effects  of  facial  features  and 
deformability  of  ocular  tissues  on  the  blast  impact.  In  addition,  we  have  made 
significant  progress  towards  developing  a  constitutive  model  for  the  stress  response 
of  the  cornea  and  sclera  that  includes  the  anisotropic  effects  of  the  collagen  fiber 
structure.  Moreover,  we  have  begun  developing  analysis  methods  for  the  dynamical 
inflation  experiments.  These  accomplishments  are  important  steps  in  examining  all 
three  hypotheses. 

2.0  BODY 

2.1  Computational  model  to  investigate  effect  of  facial  features  and  orbital 
tissue  properties  on  blast  wave  loading  on  the  eye. 

This  section  presents  a  computational  model  of  the  human  eye  that  takes  into 
account  the  propagation  of  the  wave  and  the  associated  deformation  of  the  soft 
tissue  of  the  eye  and  orbit.  The  objective  of  this  study  was  to  evaluate  the  influence 


of  facial  features  on  blast  wave  loading  on  the  human  eye,  to  simulate  deformation 
of  the  globe  arising  from  realistic  pressure  loadings,  and  to  characterize  the  time- 
history  of  the  intraocular  pressure  (IOP),  internal  stresses,  and  deformations  caused 
by  the  propagating  blast  wave.  The  result  of  this  study  has  been  submitted  for 
publication  and  the  submitted  manuscript  is  included  in  the  Appendix  of  this  report. 

The  model  involves  three-dimensional  geometries,  moving  structure 
boundaries  within  a  fluid  domain,  and  large  flow-induced  deformations  of  the 
structure.  The  flow  field  around  the  structure  is  compressible  and  pressure  fields  on 
the  structure  are  non-uniform  and  highly  transient.  The  modeling  of  the  structure 
involves  geometric  and  material  nonlinearity.  The  coupled  fluid-structure 
interaction  (FSI)  solver  uses  the  following  methods:  (1)  finite  difference 
compressible  flow  solver  for  the  propagation  of  the  blast  wave,  (2)  finite  element 
elastodynamic  solver  with  finite  deformation  of  the  globe,  and  (3)  a  sharp-interface 
immersed  boundary  method  for  fluid-structure  interaction.  The  flow  computations 
were  performed  on  an  Eulerian  grid  while  the  immersed  structural  surfaces  were 
tracked  in  a  Lagrangian  framework.  The  detail  development  of  each  method  is 
provided  for  in  the  Appendix. 

2.1.1  Finite  Element  Models  of  the  Human  Face  and  Eye 

To  characterize  pressure  loading  on  the  eye,  we  first  modeled  the  interaction 
of  the  blast  wave  with  a  rigid  (non-deformable)  finite  element  model  of  the  human 
head.  The  original  head  model,  including  the  skin  and  skeletal  structure  among 
other  anatomical  components,  was  developed  by  Zygote  Media  Group,  Inc.©[l].  This 
3D  triangular  surface  mesh  was  modified  at  the  US  Army  Research  Laboratory 
(ARL)  using  a  suite  of  advanced  mesh  processing  software,  and  quadratic  edge 
decimation  to  form  a  smooth,  manifold  surface.  Altair  HyperMesh©,  a  finite  element 
pre-processor,  was  used  to  convert  the  surface  to  a  volumetric  finite  element 
tetrahedral  mesh  [2].  In  this  development,  the  skin  was  meshed  through  to  the 
skull;  therefore,  nodes  are  shared  at  the  interface.  The  skin  mesh,  referred  to  as  the 
"skin  model"  in  the  following  developments,  is  shown  in  Figure  1A.  The  skin  model 
was  scaled  to  the  dimensions  of  a  specific  21  year  old  male  (similar  to  50  percentile 


male)  to  match  the  orbital  width  and  orbital  height  as  given  in  Weaver  et  al.  [3].  The 
interpupillary  distance  of  this  subject  was  measured  as  68.3  mm.  The  dimensions  of 
the  skin  model  are  shown  in  Figure  1B,C. 


055 


055 


Figure  1:  FE  mesh  of  the  skin  model 


To  characterize  the  deformation  of  the  globe  resulting  from  the  blast  wave,  we 
considered  the  skull  as  a  rigid  structure  (Figure  2).  The  eyelid  and  skin  were 
neglected  in  this  model,  which  is  referred  to  as  the  "skull  model". 


Skull  Orbit  Globe 


Cross  section  of  the  orbit  (Section  AB) 


Figure  2:  Skull  model:  FE  model  of  the  skull,  orbit  and  the  globe  used  in  investigating 
the  biomechanical  response  of  the  globe. 


Top  view 


Figure  3:  Definition  of  Lateral  eye  protrusion  (LP)  and  lateral  distance  (LD). 


The  model  of  the  globe  was  constructed  to  fit  within  the  orbit  of  the  skull  as  shown 
in  Figure  2,  based  on  the  anatomical  measurements  obtained  using  CT  scans  by 
Weaver  et  al.  [3].  Two  parameters,  namely  the  lateral  eye  protrusion  (LP)  and  the 
lateral  distance  (LD),  were  used  by  Weaver  et  al.  [3]  to  characterize  the  anatomical 
measurements  of  different  subjects.  The  anterior-posterior  axis  was  rotated  in  the 
axial  plane  of  the  globe  with  maximum  eye  protrusion  until  it  was  aligned  with  the 
cornea  and  the  center  of  the  optic  canal.  A  normal  axis  to  the  anterior-posterior  axis 
was  drawn  on  the  edge  of  the  lateral  orbital  rim  ([3],  Figure  2).  Lateral  eye 
protrusion  (LP)  and  lateral  distance  (LD)  are  defined  as  the  distance  of  the  cornea  to 
the  axes  intersection  and  distance  of  axes  intersection  to  the  lateral  orbital  rim, 
respectively.  To  place  the  globe  in  the  orbit,  LP  and  LD  values  were  matched  as  in 
Weaver  et  al.  [3]  for  a  specific  21  year  old  male  (similar  to  50  percentile  male).  The 
anthropometric  parameters  are  given  in  Table  2. 


Anatomical  parameter 

Value  (mm) 

Globe  diameter 

25.0 

Orbital  width,  OW 

36.0 

Orbital  height,  OH 

29.9 

Lateral  eye  protrusion,  LP 

12.0 

Lateral  distance,  LD 

19.0 

Interpupillary  distance,  IPD 

68.3 

Table  1:  Anatomical  parameters  used  in  skin  and  skull  model 


The  finite  element  mesh  of  the  globe  was  created  in  CUBIT©,  a  geometry  and  mesh 
generation  toolkit  developed  at  Sandia  National  Laboratories  [4].  The  mesh  is 
composed  of  linear  hexahedral  elements.  The  globe  included  a  corneo-scleral  shell 
and  a  space-filling,  homogeneous  internal  solid  (hereafter  referred  to  as  "inner 
solid")  in  lieu  of  vitreous  and  aqueous  humor  (Figure  2).  The  model  did  not  include 
internal  ocular  components,  namely,  lens,  zonules,  and  ciliary  muscle.  The  diameter 
of  the  globe  was  taken  as  25  mm  [5]  and  its  cross-section  was  comprised  of  two 
spherical  geometries,  namely,  the  anterior  and  the  posterior  chamber.  The  posterior 
chamber  represented  the  sclera  with  an  inner  and  outer  radius  of  12  and  12.5  mm, 
respectively.  The  anterior  chamber  represented  the  cornea,  formed  by  a  7.8  mm 
radius  circle  whose  center  was  offset  by  5  mm  from  the  center  of  the  posterior 
chamber  [5].  The  outer  and  inner  radii  of  the  cornea  were  7.8  and  7.3  mm, 
respectively.  As  a  first  approximation,  the  thickness  of  the  corneo-scleral  shell  was 
considered  uniform  with  a  value  of  0.5  mm.  The  globe  was  embedded  in  a  cone  of  a 
space-filling,  homogeneous  solid  (hereafter  referred  as  "outer  solid”)  representing 
the  extra-ocular  tissues  and  orbital  fat  (Figure  2).  The  stiffer  extra-ocular  eye 
muscles  and  internal  components  such  as  the  lens,  retina,  and  optic  nerve  head  were 
not  represented  explicitly  in  this  model.  The  attachment  point  of  the  orbital  fat  on 
the  globe  was  determined  from  the  measurements  for  a  healthy  25-year-old  female, 
reported  by  Schutte  et  al.  [6]. 

2.1.2  Material  properties 

In  our  model,  the  bulk  modulus  (K)  and  density  ( p  )  of  the  internal  solid 
were  based  on  the  respective  average  values  of  the  aqueous  and  vitreous  humor. 

The  densities  of  the  aqueous  and  vitreous  humor  were  reported  as  1003  and  1009 
kg  m-3  respectively  [7],  and  the  density  of  the  internal  solid  was  taken  as  the 
average  value  (1006  kg  m-3).  The  speed  of  sound  was  reported  as  1481-1525  m/s 
in  aqueous  at  25.5°C  and  as  1523-1532  m/s  in  vitreous  at  37°C  [7].  We  took  the 
average  speed  of  sound  in  aqueous  and  vitreous,  1503  m/s,  in  the  internal  solid  [5] 


to  calculate  the  bulk  modulus  as  2.272  x  103  MPa.  The  shear  modulus  of  the  internal 
solid  was  taken  as  7  Pa,  from  the  quasi-steady  measurements  of  the  shear  modulus 
of  bovine  vitreous  in  Ref.  [8].  The  density  of  the  corneo-scleral  shell  was  taken  as 
1400  kg  m-3  [5].  The  bulk  modulus  of  the  corneo-scleral  shell  was  taken  as  3.5706x 
103  MPa,  which  is  based  on  the  measured  speed  of  sound  in  the  sclera,  1597  m/s 
[9].  The  dynamic  shear  modulus  (G)  of  the  corneo-scleral  shell  was  taken  as  1  MPa 
from  the  low-strain  rate  measurements  of  Bisplinghoff  et  al.  [10]  for  the  sclera.  The 
low-strain  rate  values  were  chosen  to  be  consistent  with  the  available  properties  for 
the  vitreous  and  extraocular  tissues.  The  effect  of  the  scleral  shear  modulus  on  the 
deformation  of  the  globe  will  be  examined  as  a  parameter  study  in  Section  3.2.1. 

The  material  properties  of  the  outer  solid  were  taken  as  the  properties  reported  for 
the  orbital  fat.  The  density  and  bulk  modulus  of  the  outer  solid  were  assumed  to  be 
equal  to  that  of  water  [11-13]  and  taken  as  1000  kg/m3  and  2.202  x  103  MPa 
(assuming  speed  of  sound  in  water  as  1484  m/s),  respectively.  The  shear  modulus 
of  the  outer  solid  was  taken  as  500  kPa  from  the  measurements  in  Ref.  [14].  The 
properties  for  the  different  components  are  summarized  in  Table  2. 


Material 

P 

(kg  in3) 

K 

(MPa) 

G 

(Pa) 

Innersolid 

1006 

2.272  x  103 

7 

Corneo-scleral  shell 

1400 

3.571x  103 

lx  106 

Outer  solid 

1000 

2.202  x  103 

5x  105 

Table  2:  Material  parameters  used  in  the  model 


2.1.3  Modeling  the  TNT  blast  wave 

The  explosion  of  TNT  produces  a  mixture  of  gases  in  the  form  of  a  "fireball"  with 
high  temperature  and  velocities  [15].  McNesby  et  al.  [15]  measured  the  shock  front 
velocity  and  temperature  for  an  unconfined  explosion  of  2  kg  TNT  charge.  We 
applied  the  measured  velocity  and  temperature  of  the  shock  front  from  McNesby  et 
al.  [15]  as  an  initial  condition  for  the  simulations  at  a  radial  distance  of  1.8  m  from 
the  center  of  the  explosion.  The  pressure  of  the  shock-front  for  unconfined  TNT 
explosions  was  well-documented  by  Kingery  and  Bulmash  [16].  This  data  was  later 
revisited  by  Swaidak  [17]  and  is  available  in  the  form  of  scaled  distance  and  mass  of 


the  TNT.  It  has  been  widely  accepted  by  researchers  and  the  fitted  curves  have  been 
implemented  in  the  ConWep  code  of  ARL  and  LS-DYNA.  We  applied  the  pressure 
measured  for  the  shock  front  generated  by  a  2  kg  TNT  explosion  at  1.8  m  from  the 
center  of  the  explosion  [16, 17].  The  Gaussian  profiles  of  the  initial  pressure  and 
velocity  are  shown  in  Figure  4A  and  the  initial  temperature  for  the  shock  wave  is 
shown  in  Figure  4B.  The  initial  condition  elsewhere  in  the  computational  domain 
was  set  to  zero  for  velocity,  and  ambient  value  for  the  pressure  and  temperature. 


Figure  4:  Initial  conditions  for  the  pressure,  velocity  and  temperature. 


2.1.4  Results 

We  considered  a  2  kg  TNT  explosion  in  front  of  the  face  at  a  distance  of 
2.5  m  (Figure  5).  These  parameters  are  based  on  the  conditions  of  field  tests 
for  blast  exposure,  conducted  by  the  US  Army  Research  program  [15]. 


Figure  5:  Simulation  set  up 


To  simulate  the  effects  of  facial  features  on  blast  wave  reflections  around  the  eye, 
the  skin  and  all  facial  features  including  the  eye  were  treated  as  rigid  structures.  The 
flow  fields  in  the  sagittal  plane  obtained  by  the  fluid-structure  interaction 
simulations  are  shown  in  Figure  6A.  There  were  extensive  reflections  of  the  blast 
wave  by  the  nasal  and  brow  ridges.  Together,  they  formed  a  reflector  that  focused 
the  incident  blast  wave  on  the  eye  (Figure  6A).  This  effect  allows  a  planar  blast 
wave  with  a  pressure  p  just  before  impacting  the  face  to  apply  a  peak  pressure  of  4p 
(1.4  MPa)  on  the  eye  during  the  impact.  These  reflections  also  produced  an 
asymmetric  pressure  loading  on  the  eye  (Figure  6B)  that  was  higher  nearer  to  the 
nasal  bone  than  the  zygomatic  bones.  Figure  7  shows  the  time-variation  of  the 
pressure  on  the  numerical  probes  A,  B  and  C  shown  in  Figure  6B.  The  peak  pressure 
at  probe  B  was  around  two  times  higher  than  the  pressures  at  probes  A  or  C,  which 
confirms  the  asymmetric  nature  of  pressure  loading.  The  time-scale  in  Figure  7 
further  illustrates  the  highly  transient  and  asymmetric  nature  of  the  pressure 
loading  developed  from  blast  wave  interactions.  The  duration  of  pressure  loading 
on  the  face  was  less  than  1  ms  (Figure  7),  which  is  300  times  faster  than  the  blink  of 
the  eye  (~  300  ms)  [18], 


Figure  6:  (A)  Flow  field  and  pressure  contours  are  shown  in  sagittal  plane.  Zoomed 
in  view  of  the  brow  area  is  shown  to  illustrate  the  blast-wave  reflections  (B) 


Pressure  contours  on  the  face  at  the  instance  of  the  maximum  pressure  loading 
during  the  impact  the  blast  wave. 


Figure  7:  Comparison  among  time-varying  pressures  at  probes  A,  B  and  C  shown  in 
the  Fig  IB. 

The  influence  of  the  location  of  the  explosion  was  evaluated  for  a  2  kg  TNT 
explosion  on  the  ground  at  a  distance  of  2.5  m  (Figure  8).  Simulations  of  the  blast  on 
the  ground  suggested  a  similar  asymmetric  pressure  loading  on  the  eye  as  described 
in  previous  section.  The  peak  pressure  decreased  by  40%  for  the  explosion  on  the 
ground  compared  to  that  in  front  of  the  face. 


explosion 
of  2  kg  TNT 

Figure  8:  Initial  condition  of  the  simulation  for  influence  of  blast  location. 


n  sagittal  plane. 


Figure  9:  Flow  fi 


The  flow  field  shown  in  Figure  9  illustrates  that  the  nasal  and  brow  ridges  do 
not  act  as  an  effective  reflector  of  the  incident  wave  in  this  case. 

The  skull  model  described  in  section  2.2.1  was  used  to  evaluate  the  deformation 
response  of  the  eye  to  the  incident  blast  wave.  All  deformable  solids,  namely, 
internal  solid,  corneo-scleral  shell  and  outer  solid  were  treated  as  neo-Hookean 
quasi-incompressible  solids  in  the  structure  model.  The  calculated  peak  pressure 
loading  on  the  skull  model  was  30%  lower  as  compared  to  the  skin  model.  In  the 
skull  model,  removing  the  skin  layer  decreased  the  brow  protrusion  and  the  lower 
cartilaginous  part  of  the  brow,  which  decreased  the  wave  reflections  on  the  eye. 
Hence,  the  peak  pressure  was  lower  than  for  the  skin  model.  The  skull  model  also 
exhibited  asymmetric  pressure  loading  on  the  eye.  We  examined  the  question  of 
whether  the  pressure  loading  on  the  eye  changes  due  to  globe  deformation,  and 
determined  that  globe  deformation  had  a  negligible  effect  (less  than  1%)  on  the 
loading  for  these  cases. 

The  time-varying  maximum  principal  stress  (si)  of  the  corneo-scleral  shell 
and  intra-ocular  pressure  (IOP)  are  plotted  in  Figure  10A  and  si  as  well  as  IOP 
increased  overall  with  time  in  addition  to  their  periodic  behavior.  The  overall 
increase  was  due  to  the  time-varying  pressure  loading  boundary  condition  and  the 
periodic  behavior  of  the  stress  and  IOP  can  be  explained  by  the  wave  propagation 
inside  the  orbit  as  follows.  The  blast  wave  travels  inside  the  globe  and  the  outer 
solid  and  becomes  scattered  when  contacting  the  bony  wall  of  the  orbit.  This 


induces  a  periodic  behavior  of  the  intra-ocular  pressure.  The  pressure  contours  in 
the  transverse  plane  are  shown  in  Figure  10B.  The  pressure  was  calculated  as  the 
mean  of  the  three  components  of  the  principal  stress.  In  Figure  10B,  the  highest 
pressure  can  be  noticed  near  the  orbital  wall  due  to  the  wave  scattering  by  the 
orbital  wall.  We  note  that  the  pressure  is  highest  in  the  posterior  region  inside  the 
inner  solid.  The  grid  points  (red  dots)  where  the  IOP  exceeded  0.35  MPa  at  1.52  ms 
are  shown  in  the  inset  of  Figure  10A.  The  contours  of  the  von  Mises  stresses  in  the 
transverse  plane  are  plotted  in  Figure  10C,  in  which  the  largest  von  Mises  stresses 
appeared  in  the  sclera  near  the  limbus.  Due  to  asymmetric  pressure  loading  the 
globe  distorted  relative  to  the  orbital  tissue  and  resulted  in  significant  distortion  of 
the  sclera  (Figure  10C). 


Figure  10:  (a)  Time-varying  corneoscleral  principal  stress  (si)  and  maximum 
intraocular  pressure  (IOP)  (b)  Maximum  pressure  in  the  transverse  plane  of  the 


orbit.  Distortion  of  the  globe  due  to  the  blast  wave  is  shown  by  the  deformed  mesh, 
(c)  von  Mises  stress  shown  in  the  transverse  plane,  which  is  the  largest  in  the  sclera 
near  limbus. 

2.1.5  Discussion:  implictions  for  eye  injury 

As  shown  in  Figure  15A  for  the  baseline  case,  the  maximum  IOP  inside  the  globe 
reached  0.42  MPa,  corresponding  to  around  3000  mm  Hg,  which  is  two  orders  of 
magnitude  larger  than  the  15  mm  Hg  physiologic  IOP  for  a  healthy  eye.  The  von 
Mises  stresses  were  calculated  largest  in  the  sclera  wall  at  the  site  of  muscle 
attachments,  which  indicates  significant  muscle  tearing  or  rupture  inside  the  eye. 
We  assessed  the  risk  of  eye  injury  using  the  functions  published  by  Duma  and 
Kennedy  [19].  These  injury  risk  functions  were  reported  using  measurements  for 
blunt  impact  on  the  eye  by  a  projectile.  Recently,  Alphonse  et  al.  [20]  calculated  the 
injury  risk  for  an  increased  IOP  using  correlation  between  the  projectile  normalized 
energy  and  the  IOP  by  Duma  et  al.  [21].  The  calculated  percentage  of  injury  risks 
corresponding  to  the  maximum  IOP  value  0.42  MPa  for  our  simulation  are  around 
98%,  14%,  0%,  0%,  0.02%  respectively,  for  corneal  abrasion,  hyphema,  lens 
dislocation,  retinal  damage  and  globe  rupture.  We  summarize  recent  relevant 
clinical  and  experimental  studies  of  primary  ocular  blast  injury  to  provide  context 
for  our  numerical  data  and  injury  risk  analysis.  Cockerham  et  al.  [22]  studied 
combat  veteran  in-patients  (n  =  46)  with  documented  blast-induced  traumatic  brain 
injury  in  Iraq  and  Afghanistan.  The  closed-eye  injuries  such  as  corneal  abrasion, 
vitreous  haemorrhage,  retinal  detachment  and  optic  nerve  atrophy  were  found  in 
43%  of  the  patients  and  it  was  noticed  that  patients  outside  the  military  vehicles 
had  higher  levels  of  ocular  injury  than  did  those  situated  within  military  vehicles 
[22].  These  authors  attributed  this  to  closer  proximity  to  the  blast  source,  with 
attendant  high  blast  overpressure,  and  the  protective  effect  of  armoured  vehicles 
against  blast  waves  and  fragmentation.  Hines-Beard  et  al.  [23]  reported  primary 
ocular  blast  injury  on  mice,  which  were  exposed  to  a  maximum  of  0.21  MPa  blast 
pressure.  In  these  experiments,  increased  ocular  damage  was  reported  with 
increased  blast  pressures  and  the  following  closed  eye  injuries  were  reported: 
corneal  edema,  corneal  abrasions,  and  optic  nerve  avulsion. 


2.2  Dynamic  Inflation  Experiments 


We  have  designed  and  built  a  shock  tube  setup  coupled  with  full  field 
displacement  measurements  to  characterize  the  mechanical  behavior  of  the  cornea 
and  sclera  at  high  pressure  rates  representative  of  primary  blast  conditions.  The 
idea  behind  the  dynamic  inflation  experiments  is  to  record  the  dynamic  pressure 
and  full  field  displacement  history  at  various  levels  of  pressure.  We  previously 
described  the  design  and  validation  of  the  shock  tube  in  the  FY11  Annual  Report.  A 
schematic  of  the  shock  tube  setup  for  dynamic  inflation  experiments  is  shown  in 
Figure  11.  The  shock  tube,  a  specimen  holding  fixture  and  2  high-speed  imaging 


cameras  are  the  components  of  the  setup. 
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Figure  11:  Schematic  of  a  shock  tube  setup  for  dynamic  inflation  experiment. 

In  FY12,  work  has  focused  on  understanding  the  effect  of  inertia  on  the 
evolution  of  the  stress  state  of  the  membrane  specimen  subjected  to  dynamical 
inflation.  This  is  necessary  to  determine  the  stress  response  of  the  corneal  and 
sclera  specimens  under  dynamic  inflation.  We  performed  a  series  of  tests  on  an 
initially  flat  elastomeric  membrane.  The  membrane  was  composed  of  Sylgard-184, 
which  is  a  silicone  elastomer  (trademark  of  Dow  Corning)  with  well  known 
mechanical  properties.  Figure  12  plots  the  out  of  plane  displacement  of  the  top 
surface  for  a  cross-section  of  the  central  region  of  the  membrane  measured  by  DIC. 
Also  plotted  are  ellipse  fits  of  the  deformed  surface  that  were  used  to  calculate  the 
curvature  of  the  deformed  membrane  at  the  apex.  The  membrane  appears  to 


deform  rigidly  with  negligible  curvature.  The  curvature  develops  at  the  edge  of  the 
specimen  and  travels  radially  inwards  towards  the  apex.  Figure  13A  shows  a  plot  of 
the  principal  in-plane  Green-Lagrange  strain  at  the  apex  of  the  specimen.  The  strain 
field  was  calculated  from  the  gradient  of  the  displacement  field  measured  by  DIC. 
The  apex  does  not  deform  until  0.4  ms,  after  which  both  strain  components  increase 
exponentially  with  time.  Figure  13A  plots  the  pressure  as  a  function  of  apex 
displacement,  measured  by  internal  pressure  gages  near  the  specimen.  The  results 
show  a  flat  pressure-displacement  response.  In  contrast,  pressure  increases  with 
apex  displacement  for  quasistatic  inflation. 


Figure  12.  The  time  evolution  of  the  out-of-plane  displacement  response  plotted  for 
a  line  x=0  for  the  Sylgard  membrane  undergoing  dynamic  inflation. 
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Figure  13  (A)  The  principal  in-plane  strain  response  with  time  at  the  apex  of  the 
inflating  membrane,  and  (B)  the  pressure  displacement  response  for  the  apex. 


To  understand  the  pressure-displacement  response,  we  performed  a  finite 
element  simulation  of  the  dynamic  inflation  experiments.  A  ramp-hold  pressure 
loading  was  applied  to  the  inner  surface  of  the  membrane,  where  the  rise  time  and 
pressure  level  was  taken  from  pressure  measurements  of  the  shock  tube 
experiments.  Figure  14  shows  the  deformed  configuration  of  the  membrane  shortly 
after  the  pressure  loading  reached  a  constant  value.  The  simulation  shows  that 
dynamic  pressurization  sets  up  a  bending  wave  that  travels  radially  from  the 
clamped  boundaries  to  the  apex.  The  central  region  of  the  membrane  displaces 
vertically  and  remains  flat  until  the  bending  wave  passes  through  the  region  setting 
up  the  equilibrium  curvature.  In  FY13,  we  will  continue  the  computational  study  of 
the  response  of  a  flat  membrane  and  ellipsoidal  dome,  geometry  representative  of 
the  cornea  and  sclera,  to  dynamic  inflation.  The  results  of  the  studies  will  guide  the 
development  of  a  stress  analysis  for  the  experiments.  In  addition,  we  will  further 
develop  the  dynamic  inflation  experiments  to  characterize  the  properties  of  porcine 
cornea  and  scleral  cups.  This  task  will  involve  development  of  speckling  methods 
that  are  robust  enough  to  withstand  the  internal  blast  loading  and  provides 
sufficient  resolution  of  the  DIC  measurements  for  local  anisotropy. 
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Figure  14  Finite  element  simulation  of  an  initially  flat  membrane  subjected  to  a 
dynamic  constant  pressure  loading  on  the  internal  surface. 


3.0 


KEY  RESEARCH  OUTCOMES 


1.  Built  a  shock  tube  system  with  digital  image  correlation  (DIC)  setup  to 
perform  dynamic  inflation  experiments  of  ocular  tissues  and  blast 
experiments  on  the  globe. 

2.  Validated  the  shock-tube  design  by  comparing  theoretical  and  measured 
shock  pressures. 

3.  Measured  the  accuracy  of  the  DIC  displacement  measurements  for  both 
quasistatic  and  dynamic  test  systems. 

4.  Developed  a  method  for  strain  calculation  of  DIC  displacement 
measurements. 

5.  Performed  a  preliminary  dynamic  inflation  experiment  on  bovine  cornea  and 
porcine  sclera. 

6.  Developed  a  fluid-structure  interaction  solver  that  involves  compressible 
flow  and  large  deformation. 

7.  Developed  computational  model  of  human  eye  and  applied  the  model  to 
study  the  effect  of  facial  features  and  tissue  properties  on  determining  the 
blast  pressure  loading  on  the  eye. 

8.  Began  developing  a  microstructure-based  constitutive  model  to  describe  the 
nonlinear,  anisotropic  stress  response  of  the  tissues. 

9.  Measured  the  anisotropic  collagen  orientation  distribution  of  the  posterior 
human  sclera. 

4.0  Reportable  Outcomes 

1.  Bhardwaj  R.,  Ziegler  K,  Seo  J.  H,  Ramesh  K.T.,  Nguyen  T.D.  (2012)  “A 
Computational  Model  of  Blast  Loading  on  the  Human  Eye”,  Biomechanics  and 
Modeling  in  Mechanobiology,  submitted  October  2012. 


2.  Ziegler  K.A.,  Nguyen  T.D.  “Modeling  study  investigating  the  depth-dependent 
fiber-reinforcement  in  comeal  tissue”,  Poster  Presentation  Summer 
Bioengineering  Conference,  June  20-23th,  2012,  Fajardo  Puerto  Rico. 

3.  Bhardwaj  R.,  Ziegler  K.  A.,  Nguyen  T.  D.  “Blast  Wave  Reflections  on  the  Human 
Eye”,  Poster  Presentation:  ARVO,  Fort  Lauderdale,  FL,  May  9,  2012. 

4.  Ziegler  K.A.,  Yatnalkar  R.S.,  Ramesh  K.T.,  Nguyen  T.D.  “Modeling  Study  for 
the  Design  of  An  Innovative  Composite  Membrane  Inflation  Test”,  Poster 
Presentation:  Summer  Bioengineering  Conference,  June  22-25th,  2011, 
Farmington  Pennsylvania. 

5.  Yatnalkar  R.S.,  Ziegler  K.A.,  Ramesh  K.T.,  Nguyen  T.D.,  “Development  of  a 
shock  tube  to  study  primary  blast  effects  in  ocular  blast  injuries”,  Conference 
presentation:  Society  of  Experimental  Mechanics,  June  14th,  2011,  Uncasville, 
CT. 

4.0  CONCLUSION 

The  main  focus  over  the  second  year  of  the  project  was  to  develop  a  fluid- 
structure  interaction  solver  and  a  computational  model  of  the  human  eye  to 
examine  the  effects  of  facial  features  and  tissue  properties  on  the  blast  pressure 
loading.  The  results  show  that  the  nasal  and  brow  ridge  act  as  reflector  that  focuses 
the  pressure  loading  on  the  eyes.  This  also  generated  an  asymmetric  pressure 
loading  on  the  eye  that  leads  to  gross  distortion  of  the  globe  within  the  extraocular 
tissues.  The  deformation  of  the  globe  had  little  effect  on  the  blast  pressure  loading 
on  the  eye.  We  are  currently  developing  an  anisotropic  hyperelastic  constitutive 
model  for  the  tissues  of  the  eye  wall  to  determine  the  effects  of  collagen  structure  on 
the  dynamic  response  of  the  tissues.  This  will  further  advance  the  accuracy  of  the 
computational  model  of  the  eye.  We  continued  to  make  progress  on  the 
development  of  the  dynamic  inflation  test.  Work  in  the  second  year  focused  on 
understanding  the  evolution  of  the  stress  and  strain  fields  from  dynamic 
pressurization.  This  is  needed  to  develop  a  stress  analysis  for  the  tests,  which  is 


ultimately  needed  to  determine  mechanical  properties  of  the  tissues.  These 
accomplishments  add  to  the  scientific  knowledge  base  for  ocular  biomechanics 
under  blast  conditions  as  well  as  the  broader  knowledge  base  of  computational 
biomechanics  and  experimental  mechanics. 

For  FY13,  we  plan  to  complete  the  development  of  the  analytical  methods  for 
the  dynamic  inflation  tests.  Furthermore,  we  will  apply  the  method  to  measure  the 
dynamic  stress-strain  response  and  rupture  conditions  of  human  cornea  and  sclera. 
For  the  modeling  effort,  we  will  complete  development  of  a  stress-strain  model  for 
the  cornea  and  sclera  that  incorporates  the  anisotropic  lamellar  structure  and 
crimped  fibril  structure.  The  parameters  of  the  model  will  be  determined  from  the 
inflation  measurements.  The  constitutive  model  will  be  integrated  into  the  current 
computational  model  for  the  eye.  We  will  also  further  develop  the  computation 
model  to  include  important  intraocular  components,  such  as  the  lens,  retina,  and 
optic  nerve  head.  The  computational  model  will  be  applied  to  investigate  the 
deformation  and  stress  fields  in  the  internal  ocular  components.  Furthermore,  we 
will  apply  the  model  to  simulate  the  effect  of  protective  eye  equipment  on  blast 
wave  loading  on  the  eye.  These  simulations  will  ultimately  help  to  design  next 
generation  protective  eyewear,  which  could  drastically  mitigate  blast  injuries. 
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Abstract 

Ocular  injuries  from  blast  have  increased  in  recent  wars,  but  the  injury  mechanism 
associated  with  the  primary  blast  wave  is  unknown.  We  employ  a  three-dimensional 
fluid- structure  interaction  computational  model  to  understand  the  stresses  and 
deformations  incurred  by  the  globe  due  to  blast  overpressure.  Our  numerical  results 
demonstrate  that  the  blast  wave  reflections  off  the  facial  features  around  the  eye  increase 
the  pressure  loading  on  and  around  the  eye.  The  blast  wave  produces  asymmetric  loading 
on  the  eye,  which  causes  globe  distortion.  The  deformation  response  of  the  globe  under 
blast  loading  was  evaluated  and  regions  of  high  stresses  and  strains  inside  the  globe  were 
identified.  Our  numerical  results  show  that  the  blast  loading  results  in  globe  distortion 
and  large  von  Mises  stresses  in  the  sclera,  which  may  increase  the  risk  for  muscle  tearing 
and  rupture. 

1  Introduction 

Ocular  injuries  from  blasts  have  increased  in  recent  engagements  [1,  2],  Nearly  80%  of 
ocular  injuries  suffered  in  the  Iraq  war  were  caused  by  blasts  from  munitions  and 
improvised  explosive  devices  (IEDs)  [3,  4].  Such  injuries  are  classified  as  primary  and 
secondary  injuries  and  are  caused  respectively  by  blast  over  pressures  and  fragments 
from  debris  [5].  Several  factors  are  hypothesized  to  influence  the  severity  of  primary  eye 
injury  in  the  literature  [6]:  the  high  pressure  shock  front  and  the  subsequent  wave  of 
lower  sub-atmospheric  pressure  [7],  threshold  overpressure  [6],  and  reflection  of  the 
shock  wave  by  the  orbit  [8,  9,  10].  There  is,  however,  a  dearth  of  clinical  data  that  could 
verify  the  hypotheses  above  and  establish  the  mechanism  of  the  injury.  Measuring  and 
assessing  the  influence  of  these  factors  is  difficult  because  survivable  primary  blast 
injuries  are  likely  accompanied  by  injuries  from  fragments  and  blunt-force  trauma  and 
are  thus  more  difficult  to  distinguish  and  enumerate.  Moreover,  the  severity  of  the  blast 
injuries  and  distance  of  the  tertiary  care  facility  from  the  injury  site  means  that  often 
patients  are  unable  to  recount  the  injury  event,  and  witnesses  are  unavailable.  The  same 
limitations  hinder  clinical  studies  of  the  effectiveness  of  current  eye  armor,  developed  for 
ballistic  and  laser  protection,  in  preventing  blast  injuries. 

A  few  computational  models  have  been  developed  to  study  the  effects  of  impact 
from  fragments  of  comparable  size  to  the  eye  (blunt  impact).  The  3D  model  of  Uchio  et 
al.  [11]  and  the  axisymmetric  model  of  Stitzel  et  al.  [12]  focused  on  achieving 
physiologic  geometric  representations  of  the  cornea,  sclera,  lens,  ciliary  body,  choroid, 
retina,  aqueous  humor  and  vitreous  humor.  Neither  model  considered  the  extraorbital 
tissues  and  skull  beyond  applying  a  fixed  boundary  condition  at  the  posterior  surface. 
Both  models  idealized  the  material  behavior  of  the  ocular  tissues  as  nearly 
incompressible,  linear  elastic,  and  isotropic,  though  the  model  of  Stitzel  et  al.  [12] 
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allowed  for  large  deformations.  Recognizing  the  important  structural  function  of  the  eye- 
wall,  Uchio  et  al.  [11]  developed  a  more  representative  nonlinear  constitutive  model  for 
the  stress  response  of  the  cornea  and  sclera.  The  material  properties  and  failure  strains  of 
the  cornea  and  sclera  were  measured  from  uniaxial  strip  tests  under  quasi-static  (low  rate) 
loading  conditions.  Stitzel  et  al.  [12]  employed  a  finite-element  model  to  simulate  the 
impact  of  BB  pellets,  baseballs,  and  foam  projectiles  on  the  cornea.  The  results  showed 
the  highest  principal  stresses  occurred  near  the  limbus  for  the  impacting  foam  and  near 
the  equator  for  the  baseball,  which  agreed  with  experiments  [12].  In  more  recent  work 
[13],  the  orbit  and  extraocular  tissues  were  modeled  using  a  chamber  filled  with  a 
compliant  material  and  threshold  values  of  intraocular  pressure  (IOP)  for  globe  rupture 
and  peak  values  of  stresses  in  corneoscleral  shell  were  reported.  Weaver  et  al.  [14] 
evaluated  eye  injury  due  to  blunt  impact  for  different  orbital  anthropometries  using  a 
model  of  the  globe  developed  in  [12].  Their  numerical  results  suggested  the  eye  is  more 
protected  from  impact  with  greater  brow  protrusion,  less  eye  protrusion,  and  smaller 
orbital  apertures,  provided  that  the  aperture  is  large  enough  to  deter  contact  between  the 
orbit  and  the  eye  [14].  Computational  FE  models  have  also  been  used  to  assess  the  injury 
to  the  orbit  [15,  16]  and  the  optic  nerve  [17]  as  a  result  of  blunt  object  impact  to  the  eye. 
To  the  best  of  our  knowledge,  there  is  no  model  in  the  literature  which  simulates  the  blast 
wave  outside  the  eye  and  the  associated  deformation  of  the  eye.  To  this  end,  we  present  a 
computational  model  that  takes  into  account  the  propagation  of  the  wave  and  the 
associated  eye  deformation  (with  consideration  of  extraocular  tissues  of  the  orbit).  The 
objective  of  this  study  was  to  evaluate  the  influence  of  facial  features  on  blast  wave 
loading  on  the  human  eye,  to  simulate  deformation  of  the  globe  arising  from  realistic 
pressure  loadings,  and  to  characterize  the  time-history  of  the  intraocular  pressure  (IOP), 
internal  stresses,  and  deformations  caused  by  the  propagating  blast  wave. 

2  Methods 

The  methods  used  for  the  computational  modeling  of  blast  wave  interactions  with  the  eye 
are  presented  in  this  section.  The  modeling  involves  three-dimensional  geometries, 
moving  structure  boundaries  within  a  fluid  domain,  and  large  flow-induced  deformations 
of  the  structure.  The  flow  field  around  the  structure  is  compressible  and  pressure  fields  on 
the  structure  are  non-uniform  and  highly  transient.  The  modeling  of  the  structure  involves 
geometric  and  material  nonlinearity.  We  propose  a  coupled  fluid-structure  interaction 
(FSI)  solver  that  uses  the  following  methods: 

•  Finite  difference  compressible  flow  solver  for  the  propagation  of  the  blast  wave 

•  Finite  element  elastodynamic  solver  with  finite  deformation  of  the  globe 

•  Sharp-interface  immersed  boundary  method  for  fluid-structure  interaction 

The  flow  computations  were  performed  on  an  Eulerian  grid  while  the  immersed  structural 
surfaces  were  tracked  in  a  Lagrangian  framework.  In  the  following  sections,  we  describe 
these  three  components  of  the  modeling. 
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2.1  Flow  solver 


To  resolve  the  propagation  and  scattering  of  a  blast  (shock)  wave,  we  considered  the  full 
compressible  Navier-Stokes  equations  for  air.  The  equations  are  written  in  a  conservative 
form  as 
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dt  dxj  ’ 
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where/?,  uh  p,  and  e  are  the  density,  velocity,  pressure,  and  total  energy,  respectively,  and 
Pj  is  the  viscous  stress,  qj  is  the  heat  flux,  and  y  is  the  specific  heat  ratio  (1.4  for  air). 
Equations  (1-3)  were  spatially  discretized  by  a  sixth-order  central  compact  finite 
difference  scheme  [18]  and  integrated  in  time  using  a  four-stage  Runge-Kutta  method. 
An  eight-order  implicit  spatial  filtering  proposed  by  Gaintonde  et  al.  [19]  was  applied  at 
the  end  of  each  time  step  to  suppress  high  frequency  dispersion  errors.  In  order  to  resolve 
the  discontinuity  in  the  flow  variables  caused  by  a  shock  wave  with  the  current  non- 
dissipative  numerical  scheme,  the  artificial  diffusivity  method  proposed  by  Kawai  and 
Lele  [20]  had  been  applied.  In  this  method,  the  viscous  stress  and  heat  flux  are  written  as 
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where  u  and  k  are  the  physical  viscosity  and  thermal  diffusivity,  respectively,  while  p*  is 
the  artificial  shear  viscosity,  /?  is  the  artificial  bulk  viscosity,  and  k  is  the  artificial 
thermal  diffusivity.  On  the  non-uniform  Cartesian  grid,  these  artificial  diffusivities  were 
adaptively  and  dynamically  evaluated  by 
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whereC,,,  C„,  and  C.  are  user-specified  constants,  c  is  the  speed  of  sound,  S  is  the 
magnitude  of  the  strain  rate  tensor, 
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A Xk  is  the  grid  spacing,  and  over-bar  denotes  Gaussian  fdtering  [21].  We  used  C, =0.002, 
C,=  l  .0,  and  C=0.01  as  suggested  in  Kawai  and  Lele  [20]  and  the  fourth  derivatives  were 
computed  by  a  fourth-order  central  compact  scheme  [18].  From  equations  (7)  and  (8), 
artificial  diffusivities  are  significantly  larger  only  in  the  region  where  the  steep  gradient 
of  flow  variables  exists,  and  ensure  numerical  stability  in  that  region. 


2.2  Structural  solver 


The  displacement  vector  d(x,  t)  describes  the  motion  of  each  point  in  the  defonned  solid 
as  a  function  of  space  x  and  time  t.  The  deformation  gradient  tensor  F±  can  be  defined  in 

dd 

terms  of  the  displacement  gradient  tensor — Las  follows: 

dxk 


Fik  =  dik  + 


ddt 
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whered/A-is  the  Kronecker  delta  symbol,  defined  as  follows: 
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The  right  Cauchy  green  tensor  is  defined  in  terms  of  the  deformation  gradient  tensor  as 
follows: 


(12) 


The  invariants  of  the  right  Cauchy  green  tensor  are  defined  as  follows: 

/j  =  \  +  4  +  A,  (13) 

h  =  \  K  +  KK  + 
h  =  \  W, 

where  A;  are  eigen  values  of  the  right  Cauchy  green  tensor.  The  strain-energy  density 
function  tjj  of  a  neo-Hookean,  quasi-incompressible  solid  is  written  as  [22]: 


where  G  and  K  are  shear  and  bulk  moduli  and 
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.  The  Cauchy  stress  is  given  in 


(15) 


where  J  =  det(F)  denotes  the  volume  change  ratio.  The  governing  equations  for  the 
structure  are  the  Navier  equations  (momentum  balance  equation  in  Lagrangian  form)  and 
are  written  as: 

(16) 


where  i  and  j  range  from  1  to  3,  ps  is  the  density  of  the  structure,  <7,  is  the  displacement 
component  in  the  i  direction,  t  is  the  time,  a  is  the  Cauchy  stress  tensor  and  f  is  the  body 
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force  component  in  the  i  direction.  The  momentum  balance  equation  was  solved  by  finite 
elements  using  the  Galerkin  method  for  spatial  discretization,  which  yielded  the 
following  system  of  ordinary  differential  equations  for  the  nodal  displacement  vector  d: 


A7  dn+ 1 H"  K  dn+ 1  —  F „ 


(17) 


where  M  is  the  lumped  mass  matrix  and  K  is  the  stiffness  matrix.  The  Galerkin  method 
was  implemented  in  Tahoe-,  an  open-source,  Lagrangian,  three-dimensional,  finite- 
element  solver  [23].  The  central-difference  method  was  used  for  the  time  integration, 
which  resulted  in  an  explicitly  and  conditionally  stable  second  order  scheme  [24]: 
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The  constraints  of  the  time-step  of  the  governing  equations  for  the  fluid  and  solid  are 
different  and  are  governed  by  the  wave  speed  in  the  respective  domain.  In  the  present 
case,  the  wave  speed  inside  the  eye  is  larger  than  that  in  the  ambient  air  (see  Section 
2.2.2)  because  of  the  higher  wave  speed  in  the  water-like  fluids  (aqueous  and  vitreous 
humor)  inside  the  eye.  Thus,  the  time-step  was  chosen  to  resolve  the  longitudinal  wave 
speed  within  the  structure. 


2.2.1  FE  model  of  the  orbit  and  the  eye 

To  characterize  pressure  loading  on  the  eye,  we  first  modeled  the  interaction  of  the  blast 
wave  with  a  rigid  (non-deformable)  finite  element  model  of  the  human  head.  The 
original  head  model,  including  the  skin  and  skeletal  structure  among  other  anatomical 
components,  was  developed  by  Zygote  Media  Group,  Inc.°[25].  This  3D  triangular 
surface  mesh  was  modified  at  the  US  Army  Research  Laboratory  (ARL)  using  a  suite  of 
advanced  mesh  processing  software,  and  quadratic  edge  decimation  to  form  a  smooth, 
manifold  surface.  Altair  HypcrMcsh  ,  a  finite  element  pre-processor,  was  used  to  convert 
the  surface  to  a  volumetric  finite  element  tetrahedral  mesh  [26].  In  this  development,  the 
skin  was  meshed  through  to  the  skull;  therefore,  nodes  are  shared  at  the  interface.  The 
skin  mesh,  referred  to  as  the  “skin  model”  in  the  following  developments,  is  shown  in 
Figure  1A.  The  skin  model  was  scaled  to  the  dimensions  of  a  specific  21  year  old  male 
(similar  to  50  percentile  male)  to  match  the  orbital  width  and  orbital  height  as  given  in 
Weaver  et  al.  [27].  The  interpupillary  distance  of  this  subject  was  measured  as  68.3  mm. 
The  dimensions  of  the  skin  model  are  shown  in  Figure  1B,C. 
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Figure  1:  FE  mesh  of  the  skin  model 


To  characterize  the  deformation  of  the  globe  resulting  from  the  blast  wave,  we  considered 
the  skull  as  a  rigid  structure  (Figure  2).  The  eyelid  and  skin  were  neglected  in  this  model, 
which  is  referred  to  as  the  “skull  model”. 


Skull  Orbit  Globe 


Cross  section  of  the  orbit  (Section  AB) 


Figure  2:  Skull  model:  FE  model  of  the  skull,  orbit  and  the  globe  used  in  investigating  the  biomechanical 

response  of  the  globe. 
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Top  view  ”  Side  view 

Figure  3:  Definition  of  Lateral  eye  protrusion  (LP)  and  lateral  distance  (LD). 


The  model  of  the  globe  was  constructed  to  fit  within  the  orbit  of  the  skull  as  shown  in 
Figure  2,  based  on  the  anatomical  measurements  obtained  using  CT  scans  by  Weaver  et 
al.  [27].  Two  parameters,  namely  the  lateral  eye  protrusion  (LP)  and  the  lateral  distance 
(LD),  were  used  by  Weaver  et  al.  [27]  to  characterize  the  anatomical  measurements  of 
different  subjects.  The  anterior-posterior  axis  was  rotated  in  the  axial  plane  of  the  globe 
with  maximum  eye  protrusion  until  it  was  aligned  with  the  cornea  and  the  center  of  the 
optic  canal.  A  normal  axis  to  the  anterior-posterior  axis  was  drawn  on  the  edge  of  the 
lateral  orbital  rim  ([27],  Figure  2).  Lateral  eye  protrusion  (LP)  and  lateral  distance  (LD) 
are  defined  as  the  distance  of  the  cornea  to  the  axes  intersection  and  distance  of  axes 
intersection  to  the  lateral  orbital  rim,  respectively.  To  place  the  globe  in  the  orbit,  LP  and 
LD  values  were  matched  as  in  Weaver  et  al.  [27]  for  a  specific  21  year  old  male  (similar 
to  50  percentile  male).  The  anthropometric  parameters  are  given  in  Table  2. 


Anatomical  parameter 

Value  (mm) 

Globe  diameter 

25.0 

Orbital  width,  OW 

36.0 

Orbital  height,  OH 

29.9 

Lateral  eye  protrusion,  LP 

12.0 

Lateral  distance,  LD 

19.0 

Interpupillary  distance,  IPD 

68.3 

Table  1 :  Anatomical  parameters  used  in  skin  and  skull  model 

The  finite  element  mesh  of  the  globe  was  created  in  CUBIT®,  a  geometry  and  mesh 
generation  toolkit  developed  at  Sandia  National  Laboratories  [28].  The  mesh  is  composed 
of  linear  hexahedral  elements.  The  globe  included  a  corneo-scleral  shell  and  a  space¬ 
filling,  homogeneous  internal  solid  (hereafter  referred  to  as  “inner  solid”)  in  lieu  of 
vitreous  and  aqueous  humor  (Figure  2).  The  model  did  not  include  internal  ocular 
components,  namely,  lens,  zonules,  and  ciliary  muscle.  The  diameter  of  the  globe  was 
taken  as  25  mm  [12]  and  its  cross-section  was  comprised  of  two  spherical  geometries, 
namely,  the  anterior  and  the  posterior  chamber.  The  posterior  chamber  represented  the 
sclera  with  an  inner  and  outer  radius  of  12  and  12.5  mm,  respectively.  The  anterior 
chamber  represented  the  cornea,  formed  by  a  7.8  mm  radius  circle  whose  center  was 
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offset  by  5  mm  from  the  center  of  the  posterior  chamber  [12].  The  outer  and  inner  radii  of 
the  cornea  were  7.8  and  7.3  mm,  respectively.  As  a  first  approximation,  the  thickness  of 
the  corneo-scleral  shell  was  considered  unifonn  with  a  value  of  0.5  mm.  The  globe  was 
embedded  in  a  cone  of  a  space-filling,  homogeneous  solid  (hereafter  referred  as  “outer 
solid”)  representing  the  extra-ocular  tissues  and  orbital  fat  (Figure  2).  The  stiffer  extra¬ 
ocular  eye  muscles  and  internal  components  such  as  the  lens,  retina,  and  optic  nerve  head 
were  not  represented  explicitly  in  this  model.  The  attachment  point  of  the  orbital  fat  on 
the  globe  was  determined  from  the  measurements  for  a  healthy  25-year-old  female, 
reported  by  Schutte  et  al.  [29]. 

2.2.2  Material  properties 

In  our  model,  the  bulk  modulus  ( K)  and  density  (p)  of  the  internal  solid  were  based  on  the 
respective  average  values  of  the  aqueous  and  vitreous  humor.  The  densities  of  the 
aqueous  and  vitreous  humor  were  reported  as  1003  and  1009  kg  m"3  respectively  [30], 
and  the  density  of  the  internal  solid  was  taken  as  the  average  value  (1006  kg  m'3).  The 
speed  of  sound  was  reported  as  1481-1525  m/s  in  aqueous  at  25.5°C  and  as  1523-1532 
m/sin  vitreous  at  37°C  [30].  We  took  the  average  speed  of  sound  in  aqueous  and  vitreous, 
1503  m/s,  in  the  internal  solid  [12]  to  calculate  the  bulk  modulus  as  2.272  x  103  MPa. 
The  shear  modulus  of  the  internal  solid  was  taken  as  7  Pa,  from  the  quasi-steady 
measurements  of  the  shear  modulus  of  bovine  vitreous  in  Ref.  [31].  The  density  of  the 
comeo-scleral  shell  was  taken  as  1400  kg  m'3  [12].  The  bulk  modulus  of  the  corneo¬ 
scleral  shell  was  taken  as  3.5706x  103  MPa,  which  is  based  on  the  measured  speed  of 
sound  in  the  sclera,  1597  m/s  [32].  The  dynamic  shear  modulus  (G)  of  the  comeo-scleral 
shell  was  taken  as  1  MPa  from  the  low-strain  rate  measurements  of  Bisplinghoff  et  al. 
[33]  for  the  sclera.  The  low-strain  rate  values  were  chosen  to  be  consistent  with  the 
available  properties  for  the  vitreous  and  extraocular  tissues.  The  effect  of  the  scleral  shear 
modulus  on  the  deformation  of  the  globe  will  be  examined  as  a  parameter  study  in 
Section  3.2.1.  The  material  properties  of  the  outer  solid  were  taken  as  the  properties 
reported  for  the  orbital  fat.  The  density  and  bulk  modulus  of  the  outer  solid  were  assumed 
to  be  equal  to  that  of  water  [34,  35,  36]  and  taken  as  1000  kg/m3  and  2.202  x  103  MPa 
(assuming  speed  of  sound  in  water  as  1484  m/s),  respectively.  The  shear  modulus  of  the 
outer  solid  was  taken  as  500  kPa  from  the  measurements  in  Ref.  [37].  The  properties  for 
the  different  components  are  summarized  in  Table  2. 


Material 

P 

(kg  in3) 

K 

(MPa) 

G 

(Pa) 

Innersolid 

1006 

2.272  x  103 

7 

Corneo-scleral  shell 

1400 

3.571x  103 

lx  106 

Outer  solid 

1000 

2.202  x  103 

5x  105 

Table  2:  Material  parameters  used  in  the  model 

2.3  Flu  id -structure  interaction  coupling 

A  partitioned  approach  was  used  to  couple  the  flow  and  the  structure  solvers  [38].  In  this 
approach,  flow  and  structure  solvers  are  coupled  such  that  they  exchange  data  at  each 


8 


time  step  (Figure  4A).  In  general,  there  are  two  coupling  methods  used  in  fluid  structure 
interaction  algorithms —  explicit  (or  weak,  one-way)  coupling  and  implicit  (or  strong, 
two-way)  coupling.  As  the  name  suggests,  explicit  and  implicit  coupling  integrate  the 
governing  equations  of  the  flow  and  the  structure  domain  explicitly  and  implicitly  in 
time.  Explicit  coupling  is  computationally  inexpensive  and  may  be  subject  to  stability 
constraints,  which  depends  on  the  structure-fluid  density  ratio  {pjpt)  [39].  On  the  other 
hand,  implicit  coupling  is  robust,  computationally  expensive  and  does  not  introduce 
stability  constraints.  Explicit  coupling  is  a  good  candidate  in  cases  where  ps/pf  is  large, 
for  example  air-tissue  interaction  during  phonation  of  vocal  folds  in  the  larynx,  while  an 
implicit  scheme  is  needed  for  low  values  of  ps/pf,  for  example  blood-tissue  interaction  in 
cardiovascular  flows.  In  the  latter  case,  the  structure  will  respond  strongly  even  with 
small  perturbations  from  the  fluid  and  vice  versa.  In  the  present  paper,  pjpf  ~  800  and 
explicit  coupling  is  used  for  the  simulations. 


A 


B 


C 


Pressure  on  the  surface 


Figure  4:  (A)  Partitioned  approach  (B)  Data  exchange  between  flow  and  solid  solver  (C)  Algorithm  of  FSI 

solver 


In  explicit  coupling,  the  flow  solution  is  marched  by  one  time  step  with  the  current 
deformed  shape  of  the  structure  and  the  velocities  of  the  fluid- structure  interface  act  as 
the  boundary  conditions  in  the  flow  solver  (Figure  4B).  This  boundary  condition 
represents  continuity  of  the  velocity  at  the  interface  (no  slip  on  the  solid  surface): 


(19) 

The  pressure  loading  on  the  structure  surface  exposed  to  the  fluid  domain  is  calculated  at 
the  current  location  on  the  structure  using  the  interpolated  normal  fluid  pressure  at  the 
boundary  intercept  points  via  a  trilinear  interpolation  (bilinear  interpolation  for  2D)  as 
described  in  Ref.  [40].  This  boundary  condition  represents  continuity  of  the  traction  at  the 
solid-fluid  interface: 


^T/.fluid^/  ®  ij  .structure  ^  j 


(20) 
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where  n}  is  the  local  surface  normal  pointing  outward  from  the  surface.  The  structure 
solver  is  marched  by  one  time  step  with  the  updated  fluid  dynamic  forces  (Figure  4c). 

2.3.1  Immersed  boundary  method 

The  compressible  Navier-Stokes  equations  for  the  fluid  flow  with  complex  structure 
boundaries  inside  the  fluid  domain  were  solved  using  the  sharp-interface  immersed 
boundary  method  of  Mittal  et  al.  [40].  In  this  method,  the  surface  of  the  structure  and  the 
fluid  domain  are  represented  by  an  unstructured  surface  mesh  and  a  Cartesian  grid, 
respectively.  The  surface  mesh  of  the  structure  and  Cartesian  grid  of  the  fluid  domain 
consists  of  triangular  elements  and  cells  (cube  or  cuboids),  respectively.  The  surface 
mesh  is  “immersed”  inside  the  fluid  grid.  At  the  pre-processing  stage  before  integrating 
governing  equations,  the  cells  of  the  fluid  domain  were  marked  according  to  their 
location  with  respect  to  the  surface  mesh.  The  cells  whose  centers  were  located  inside  the 
surface  mesh  were  identified  and  tagged  as  “body”  cells  and  the  other  points  outside  the 
surface  mesh  were  “fluid”  cells  as  shown  in  Figure  5.  Note  that  only  cell  centers  are 
shown  in  Figure  5.  Any  body-cell  which  has  at  least  one  fluid-cell  neighbor  was  tagged 
as  a  ‘ghost-cell’  (Figure  5).  The  center  of  this  ghost  cell  is  referred  to  as  a  "ghost  point". 
The  boundary  condition  on  the  surface  mesh  was  imposed  by  specifying  an  appropriate 
value  at  the  ghost-point.  A  “normal  probe”  was  extended  from  the  ghost  point  to  intersect 
the  surface  mesh.  This  intersection  point  is  referred  to  as  the  “body  intercept”.  The  probe 
was  further  extended  into  the  fluid  to  the  “image  point”  such  that  the  body-intercept  lies 
midway  between  the  image  and  ghost  points.  A  linear  interpolation  was  used  along  the 
normal  probe  to  compute  the  value  at  the  ghost  point  based  on  the  boundary-intercept 
value  and  the  value  estimated  at  the  image-point.  The  value  at  the  image-point  itself  was 
computed  through  a  tri-linear  interpolation  from  the  surrounding  fluid  nodes.  This 
procedure  leads  to  a  nominally  second-order  accurate  specification  of  the  boundary 
condition  on  the  surface  mesh  [40]. 


Figure  5:  Schematic  of  the  immersed  boundary  method 

2.3.2  Initial  conditions  associated  with  the  blast 

The  explosion  of  TNT  produces  a  mixture  of  gases  in  the  fonn  of  a  "fireball"  with  high 
temperature  and  velocities  [41].  McNesby  et  al.  [41]  measured  the  shock  front  velocity 
and  temperature  for  an  unconfined  explosion  of  2  kg  TNT  charge.  In  this  experiment,  the 
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shock  wave  separated  from  the  fireball  after  1  ms  of  the  explosion.  The  radius  of  the 
fireball  was  measured  1.8  m  at  this  instance.  We  applied  the  measured  velocity  and 
temperature  of  the  shock  front  from  McNesby  et  al.  [41]  as  an  initial  condition  for  the 
simulations  at  a  radial  distance  of  1.8  m  from  the  center  of  the  explosion.  Figure  6  shows 
the  shock  front  velocity  as  a  function  of  the  distance  for  2  kg  TNT.  The  pressure  of  the 
shock-front  for  unconfined  TNT  explosions  was  well-documented  by  Kingery  and 
Bulmash  [42].  This  data  was  later  revisited  by  Swaidak  [43]  and  is  available  in  the  form 
of  scaled  distance  and  mass  of  the  TNT.  It  has  been  widely  accepted  by  researchers  and 
the  fitted  curves  have  been  implemented  in  the  ConWep  code  of  ARL  and  LS-DYNA. 
We  applied  the  pressure  measured  for  the  shock  front  generated  by  a  2  kg  TNT  explosion 
at  1.8  m  from  the  center  of  the  explosion  [42,  43].  The  Gaussian  profiles  of  the  initial 
pressure  and  velocity  are  shown  in  Figure  7A  and  the  initial  temperature  for  the  shock 
wave  is  shown  in  Figure  7B.  The  initial  condition  elsewhere  in  the  computational  domain 
was  set  to  zero  for  velocity,  and  ambient  value  for  the  pressure  and  temperature. 


Figure  6:  Shock  front  velocity  as  a  function  of  distance.  Points  represent  measurements  by  McNesbyet  al. 

[41]. 
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3  Results  and  Discussion 


We  considered  a  2  kg  TNT  explosion  in  front  of  the  face  at  a  distance  of  2.5  m  (Figure 
8).  These  parameters  are  based  on  the  conditions  of  field  tests  for  blast  exposure, 
conducted  by  the  US  Army  Research  program  [41].  The  fluid-structure  interaction  (FSI) 
solver  presented  in  the  previous  section  was  employed  to  simulate  the  propagation  of  the 
three-dimensional  blast  wave.  The  grid  stretching  boundary  condition  was  applied  on  all 
boundaries  of  the  computational  domain  shown  in  Figure  8.  This  section  is  organized  as 
follows.  First,  we  present  the  results  of  pressure  loading  using  the  rigid  skin  model  in 
Section  3.1.  In  Section  3.2,  we  use  the  skull  model  to  characterize  the  deformation  of  the 
globe  caused  by  blast  pressure  loading.  Next,  we  present  sensitivity  study  of  the 
computational  model  to  material  properties  (section  3.2.1)  and  anatomical  parameters 
(section  3.2.2).  Finally,  we  predict  possible  eye  injuries  and  associated  injury  risks  using 
the  calculated  intra-ocular  pressure  (IOP)  and  von-Mises  stress  in  the  sclera  (section  3.3); 
and  discuss  limitations  of  our  model  in  section  3.4. 


Figure  8:  Simulation  set  up 


3.1  Effect  of  facial  features  on  blast  wave  reflections  around  the  eye 

For  these  simulations,  the  skin  and  all  facial  features  including  the  eye  were  treated  as 
rigid  structures.  The  flow  fields  in  the  sagittal  plane  obtained  by  the  fluid-structure 
interaction  simulations  are  shown  in  Figure  9A.  There  were  extensive  reflections  of  the 
blast  wave  by  the  nasal  and  brow  ridges.  Together,  they  formed  a  reflector  that  focused 
the  incident  blast  wave  on  the  eye  (Figure  9A).  This  effect  allows  a  planar  blast  wave 
with  a  pressure  p  just  before  impacting  the  face  to  apply  a  peak  pressure  of  4 p  (1.4  MPa) 
on  the  eye  during  the  impact.  These  reflections  also  produced  an  asymmetric  pressure 
loading  on  the  eye  (Figure  9B)  that  was  higher  nearer  to  the  nasal  bone  than  the 
zygomatic  bones.  Figure  10  shows  the  time-variation  of  the  pressure  on  the  numerical 
probes  A,  B  and  C  shown  in  Figure  9B.  The  peak  pressure  at  probe  B  was  around  two 
times  higher  than  the  pressures  at  probes  A  or  C,  which  confirms  the  asymmetric  nature 
of  pressure  loading.  The  time-scale  in  Figure  10  further  illustrates  the  highly  transient  and 
asymmetric  nature  of  the  pressure  loading  developed  from  blast  wave  interactions.  The 
duration  of  pressure  loading  on  the  face  was  less  than  1  ms  (Figure  10),  which  is  300 
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times  faster  than  the  blink  of  the  eye  (~  300  ms)  [44].  Note  that  the  thermal  front 
propagation  was  slower  as  compared  to  the  shock  wave.  The  shock  wave  travelled  around 
0.7  m  in  1  ms  to  reach  the  eye  while  the  thermal  front  propagated  only  0.25  m  in  this 
time.  This  can  be  attributed  to  the  low  thermal  conductivity  of  the  air. 


Figure  9:  (a)  Flow  field  and  pressure  contours  are  shown  in  sagittal  plane.  Velocity  vectors  are  shown  at 
every  10th  grid  point.  Zoomed  in  view  of  the  brow  area  is  shown  to  illustrate  the  blast- wave  reflections  (b) 
Pressure  contours  on  the  face  at  the  instance  of  the  maximum  pressure  loading  during  the  impact  the  blast 
wave.  The  pressure  of  the  blast  front  was  0.34  MPa  just  before  the  impact. 
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Figure  10:  Comparison  among  time-varying  pressures  at  probes  A,  B  and  C  shown  in  the  Fig  IB. 

3.1.1  Influence  of  blast  location 

The  influence  of  the  location  of  the  explosion  was  evaluated  for  a  2  kg  TNT  explosion  on 
the  ground  at  a  distance  of  2.5  m  (Figure  11).  Simulations  of  the  blast  on  the  ground 
suggested  a  similar  asymmetric  pressure  loading  on  the  eye  as  described  in  previous 
section.  The  peak  pressures  on  probe  B  shown  in  Figure  9B  are  plotted  in  Figure  12.  We 
note  a  40%  decrease  in  peak  pressure  for  the  explosion  on  the  ground  as  compared  to  that 
in  front  of  the  face.  The  flow  field  shown  in  Figure  13  illustrates  that  the  nasal  and  brow 
ridges  do  not  act  as  an  effective  reflector  of  the  incident  wave  in  this  case. 


explosion 
of  2  kg  TNT 

Figure  11:  Initial  condition  of  the  simulation  for  investigating  influence  of  blast  location. 
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Figure  12:  Comparison  between  time-varying  pressures  at  probe  B  for  the  explosions  in  front  of  face  and 
on  the  ground. 


Figure  13:  Flow  field  and  pressure  contours  are  shown  in  sagittal  plane.  Velocity  vectors  are  shown  at 
every  10th  grid  point.  The  pressure  levels  are  same  as  shown  in  Figure  9B. 


3.2  Biomechanical  model  of  the  eye  orbit  (baseline  case) 

The  skull  model  described  in  section  2.2.1  was  used  to  evaluate  the  deformation  response 
of  the  eye  to  the  incident  blast  wave.  All  deformable  solids,  namely,  internal  solid, 
comeo-scleral  shell  and  outer  solid  were  treated  as  neo-Hookean  quasi-incompressible 
solids  in  the  structure  model  (Eq.  14).  The  calculated  peak  pressure  loading  on  the  skull 
model  was  30%  lower  as  compared  to  the  skin  model.  The  time  variation  of  the  pressures 
on  probes  A,  B  and  C  is  shown  in  Figure  14  and  the  locations  of  the  probes  are  shown  in 
the  inset  of  Figure  14.  The  peak  pressure  (~  0.9  MPa)  was  recorded  on  probes  A  and  B. 
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In  the  skull  model,  removing  the  skin  layer  decreased  the  brow  protrusion  and  the  lower 
cartilaginous  part  of  the  brow,  which  decreased  the  wave  reflections  on  the  eye.  Hence, 
the  peak  pressure  was  lowered  by  36%  as  compared  to  the  skin  model.  The  skull  model 
also  exhibited  asymmetric  pressure  loading  on  the  eye.  We  examined  the  question  of 
whether  the  pressure  loading  on  the  eye  changes  due  to  globe  deformation,  and 
determined  that  globe  deformation  had  a  negligible  effect  (less  than  1%)  on  the  loading 
for  these  cases. 


Figure  14:  Comparison  among  time-varying  pressures  at  probes  A,  B  and  C  shown  in  the  inset. 

The  time-varying  maximum  principal  stress  (si)  of  the  comeo-scleral  shell  and  intra¬ 
ocular  pressure  (IOP)  are  plotted  in  Figure  15A  and  s i  as  well  as  IOP  increased  overall 
with  time  in  addition  to  their  periodic  behavior.  The  overall  increase  was  due  to  the  time- 
varying  pressure  loading  boundary  condition  and  the  periodic  behavior  of  the  stress  and 
IOP  can  be  explained  by  the  wave  propagation  inside  the  orbit  as  follows.  The  blast  wave 
travels  inside  the  globe  and  the  outer  solid  and  becomes  scattered  when  contacting  the 
bony  wall  of  the  orbit.  This  induces  a  periodic  behavior  of  the  intra-ocular  pressure.  To 
verify  this  hypothesis,  we  performed  a  first-order  time  scaling  analysis  of  the  wave 
propagation  as  follows.  The  time  period  of  the  longitudinal  wave  traveling  inside  the 
orbit  is  given  by  T  ~  L/v,  where  L  is  the  distance  of  the  cornea  to  the  orbital  wall  (=  36.36 
mm)  and  v  is  the  velocity  of  the  longitudinal  wave,  given  by, 

v  =  ((K+  l.3G)/p)m  (16) 

Where  K  and  G  are  the  bulk  and  shear  moduli,  respectively  and  p  is  the  density.  We 
considered  K»G  (Table  2)  and  averaged  the  values  of  K  and  p  for  the  inner  solid,  the 
comeo-scleral  shell  and  the  outer  solid  to  calculate  v  in  the  above  equation.  The  time 
period  calculated  by  the  scaling  analysis  was  around  0.03  ms.  The  time  period  calculated 
from  the  simulation  was  around  0.07  ms  (Figure  15  A),  which  is  on  the  same  order  of  the 
value  calculated  by  the  scaling  analysis.  We  therefore  conclude  that  the  periodic  behavior 
of  the  stress  and  maximum  IOP  is  due  to  wave  scattering  inside  the  orbit. 

The  pressure  contours  in  the  transverse  plane  are  shown  in  Figure  15B.  The 
pressure  was  calculated  as  the  mean  of  the  three  components  of  the  principal  stress.  In 
Figure  15B,  the  highest  pressure  can  be  noticed  near  the  orbital  wall  due  to  the  wave 
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scattering  by  the  orbital  wall.  We  note  that  the  pressure  is  highest  in  the  posterior  region 
inside  the  inner  solid.  The  grid  points  (red  dots)  where  the  IOP  exceeded  0.35  MPa  at 
1.52  ms  are  shown  in  the  inset  of  Figure  15  A.  The  contours  of  the  von  Mises  stresses  in 
the  transverse  plane  are  plotted  in  Figure  15C,  in  which  the  largest  von  Mises  stresses 
appeared  in  the  sclera  near  the  limbus.  Due  to  asymmetric  pressure  loading  (as  explained 
in  section  3.1),  the  globe  distorted  relative  to  the  orbital  tissue  and  resulted  in  significant 
distortion  of  the  sclera  (Figure  15C). 
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Figure  15:  (A)  Time -varying  corneoscleral  principal  stress  (si)  and  maximum  intraocular  pressure  (IOP) 

(B)  Maximum  pressure  in  the  transverse  plane  of  the  orbit.  Distortion  of  the  globe  due  to  the  blast  wave  is 
shown  by  the  deformed  mesh.  (C)  von  Mises  stress  shown  in  the  transverse  plane,  which  is  the  largest  in 
the  sclera  near  limbus. 
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3.2.1  Sensitivity  of  the  model  to  material  properties 

We  conducted  a  sensitivity  study  of  material  properties  of  the  corneo-scleral  shell  and 
outer  solid  on  the  deformation  response  of  the  globe.  The  bulk  moduli  of  the  inner  solid, 
comeo-scleral  shell,  and  outer  solid  were  not  varied  and  the  values  used  in  the  simulation 
are  given  in  Table  2.  The  shear  moduli  of  the  corneo-scleral  shell  and  outer  solid  were 
varied  as  shown  in  Table  3.  The  shear  modulus  of  the  outer  solid  was  decreased  by  4.4 
times  in  case  1  and  increased  to  a  very  high  value  in  case  2.  Case  2  represents  non- 
deformable  extra-ocular  tissues  and  represents  a  wall  boundary  condition  on  thesclera.  In 
case  3,  the  shear  modulus  of  the  outer  solid  was  kept  the  same  as  in  the  baseline  case 
while  increasing  the  shear  modulus  of  the  corneo-scleral  shell  10-fold,  which  is  on  the 
order  of  the  scleral  modulus  measured  for  high  pressure  rates  [33]. 


Cases 

G^comeo-scleralshell 

(MPa) 

G0uter  solid 

(Pa) 

Baseline 

i 

500 

Case  1 

i 

113 

Case  2 

i 

1.4e9 

Case  3 

10 

500 

Table  3:  Simulation  cases  for  sensitivity  to  material  properties. 


Figure  16:  Comparison  among  four  cases  with  different  material  properties  for  time-varying  comeo-scleral 
principal  stress,  Si  (A)  and  maximum  intraocular  pressure,  IOP  (B). 

The  comparison  of  time-varying  maximum  IOP  and  the  comeo-scleral  principal  stress,  si 
for  the  four  cases  described  above  is  shown  in  Figure  16A  and  15B,  respectively.  The 
maximum  IOP  and  comeo-scleral  stress  for  the  baseline  case  and  case  lwere  almost  the 
same.  It  indicates  that  the  decreasing  shear  modulus  of  the  outer-solid  by  4.4  times  did 
not  influence  the  results. 

The  maximum  IOP  and  corneo-scleral  stress  shown  in  the  transverse  plane  for 
case  2  are  shown  in  Figure  17.  We  notice  that  the  maximum  IOP  and  corneo-scleral 
stresses  decreased  90%  and  10%  respectively  as  compared  to  the  baseline  case  (Figure 
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17).  This  can  be  attributed  to  the  very  high  shear  modulus  of  the  outer  solid  which  acted 
as  a  wall-like  boundary  condition  and  reduced  the  scattering  of  the  travelling  wave  in  the 
orbit.  This  reduced  the  normal  deformation  of  the  posterior  globe  but  increased  the 
distortional  deformation  of  the  anterior  globe.  As  a  result,  normal  stresses  (which 
occurred  mainly  in  the  posterior  globe)  were  smaller  but  von  Mises  stresses  in  the 
anterior  comeal-scleral  shell  were  larger.  This  also  explains  why  the  region  of  maximum 
IOP  in  Fig.  17A  shifted  from  the  posterior  region  (baseline  case  in  Fig  15B)  to  the 
anterior  region  of  the  globe. 

The  sensitivity  of  the  shear  modulus  of  the  corneo-scleral  shell  was  investigated 
by  increasing  it  10  times  as  in  the  baseline  case  and  von  Mises  stresses  in  the  corneo¬ 
scleral  shell  increased  by  4  times,  as  shown  in  Figure  18. 


Figure  17:  (A)  Maximum  pressure  in  the  transverse  plane  of  the  orbit  for  case  2.  Distortion  of  the  globe  due 
to  the  blast  wave  is  shown  by  the  deformed  mesh.  (B)  von  Mises  stresses  shown  in  the  transverse  plane, 
which  is  the  largest  in  the  extra-ocular  tissue. 


Figure  18:  von  Mises  stresses  shown  in  the  transverse  and  sagittal  plane  for  case  3. 


3.2.2  Sensitivity  of  the  model  to  anatomical  parameters 

In  this  section,  we  present  the  simulations  for  understanding  the  sensitivities  to 
anatomical  parameters  in  the  model.  The  value  of  LP  was  increased  in  case  4,  keeping 
LD  same  (Table  4).  The  range  of  variation  of  LP  was  based  on  the  reported  anatomical 
measurements  for  39  subjects  by  Weaver  et  al.[14]. 


Cases 

LP 

LD 

(mm) 

(mm) 

19 


Baseline 

12 

19 

Case  4 

18 

19 

Table  4:  Anatomical  parameter  (LP)  variation  selected  for  the  sensitivity  study. 

In  Figure  19,  we  noticed  around  50%  reduction  in  the  IOP  and  corneo-scleral  stress  in 
case  4  with  the  increased  lateral  protrusion  (LP).  This  can  be  explained  by  the  change  in 
spatial  variation  of  pressure  loading  on  the  eye  due  to  change  in  anatomy  of  the  eye.  A 
comparison  of  the  pressure  field  in  the  transverse  plane  for  the  two  cases  is  shown  in 
Table  4.  The  pressure  loading  was  more  concentrated  on  the  extra-ocular  tissues  in  case  4 
as  compared  to  the  baseline  case.  The  reflections  from  the  protruded  eye  with  higher 
lateral  protrusion  (LP)  altered  the  spatial  gradient  of  pressure  on  the  eye  as  shown  in 
Figure  20. 


Figure  19:  Comparison  of  time-varying  corneo-scleral  principal  stress  (A)  and  maximum  intraocular 
pressure  IOP  (B)  betweentwo  cases  with  different  lateral  protrusion  (LP). 


Figure  20:  Comparison  of  pressure  field  outside  the  eye  between  two  cases  with  different  anatomical 
parameter  (A)  Baseline,  LP  =  12  mm  (B)  Case  4,  LP  18  mm. 
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3.3  Implications  to  the  eye  injury 


As  shown  in  Figure  15A  for  the  baseline  case,  the  maximum  IOP  inside  the  globe 
reached  0.42  MPa,  corresponding  to  around  3000  mm  Hg,  which  is  two  orders  of 
magnitude  larger  than  the  15  mm  Hg  physiologic  IOP  for  a  healthy  eye.  The  von  Mises 
stresses  were  calculated  largest  in  the  sclera  wall  at  the  site  of  muscle  attachments,  which 
indicates  significant  muscle  tearing  or  rupture  inside  the  eye.  We  assessed  the  risk  of  eye 
injury  using  the  functions  published  by  Duma  and  Kennedy  [45].  These  injury  risk 
functions  were  reported  using  measurements  for  blunt  impact  on  the  eye  by  a  projectile 
and  can  be  expressed  as  [45]: 


Inj  ury_  risk 


1 

1  +  e“-hx 


100% 


where  x  is  projectile  normalized  energy  [kJ/m2],  a  is  a  dimensionless  parameter  that 
depends  on  the  injury  type  and  parameter  b  is  expressed  in  terms  of  m2/kJ.  The  values  of 
a  and  b  for  different  injury  types  were  reported  in  Ref.  [45].  Recently,  Alphonse  et  al. 
[46]  calculated  the  injury  risk  for  an  increased  IOP  using  correlation  between  the 
projectile  normalized  energy  and  the  IOP  by  Duma  et  al.  [47].  The  calculated  percentage 
of  injury  risks  corresponding  to  the  maximum  IOP  value  0.42  MPa  for  our  simulation  are 
around  98%,  14%,  0%,  0%,  0.02%  respectively,  for  comeal  abrasion,  hyphema,  lens 
dislocation,  retinal  damage  and  globe  rupture.  We  summarize  recent  relevant  clinical  and 
experimental  studies  of  primary  ocular  blast  injury  to  provide  context  for  our  numerical 
data  and  injury  risk  analysis.  Cockerham  et  al.  [48]  studied  combat  veteran  in-patients  (n 
=  46)  with  documented  blast-induced  traumatic  brain  injury  in  Iraq  and  Afghanistan.  The 
closed-eye  injuries  such  as  corneal  abrasion,  vitreous  haemorrhage,  retinal  detachment 
and  optic  nerve  atrophy  were  found  in  43%  of  the  patients  and  it  was  noticed  that  patients 
outside  the  military  vehicles  had  higher  levels  of  ocular  injury  than  did  those  situated 
within  military  vehicles  [48].  These  authors  attributed  this  to  closer  proximity  to  the  blast 
source,  with  attendant  high  blast  overpressure,  and  the  protective  effect  of  armoured 
vehicles  against  blast  waves  and  fragmentation.  Hines-Beard  et  al.  [49]  reported  primary 
ocular  blast  injury  on  mice,  which  were  exposed  to  a  maximum  of  0.21  MPa  blast 
pressure.  In  these  experiments,  increased  ocular  damage  was  reported  with  increased 
blast  pressures  and  the  following  closed  eye  injuries  were  reported:  corneal  edema, 
comeal  abrasions,  and  optic  nerve  avulsion. 


3.4  Limitations  of  the  computational  model 

The  present  work  contains  a  number  of  limitations.  The  model  for  the  eye-wall  assumed 
a  uniform  thickness,  though  the  thickness  of  the  human  cornea  varies  from  0.5 -0.8  mm 
from  apex  to  limbus  [50]  and  the  thickness  of  the  sclera  can  vary  from  0.3- 1.2  mm.  The 
sclera  is  thickest  at  the  limbus  and  posterior  region  and  thinnest  at  the  equator  and 
peripapillary  region.  The  tissues  of  the  cornea  and  sclera  were  assumed  to  be  isotropic 
and  spatially  homogeneous,  but  X-ray  diffraction  experiments  have  demonstrated 
significant  anisotropy  in  the  limbus  [51]  and  peripapillary  regions  of  the  sclera  [52].  The 
neo-Hookean  model  furthermore  could  not  capture  the  dramatic  strain  stiffening  stress 
response  for  pressures.  We  are  currently  working  on  incorporating  more  detailed 
description  of  the  spatially  varying  thickness  and  collagen-derived  anisotropy  of  the 
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cornea  and  sclera  in  the  computational  model  to  investigate  their  effects  on  the 
deformation  and  stress  response  of  the  globe  to  blast.  The  computational  model  also  does 
not  include  important  intraocular  components  such  as  the  lens,  retina,  and  optic  nerve 
head.  These  components  will  be  included  in  future  models  to  investigate  the  risk  for  such 
injuries  as  lens  dislocation,  retinal  detachment,  and  optic  nerve  avulsion.  We  will  apply 
the  more  detailed  computational  model  to  investigate  the  ability  of  current  eye  protective 
equipment  to  reduce  the  effects  of  blast  loading  on  the  globe. 


4  Summary 

A  fluid-structure  interaction  computational  model  was  presented  for  blast  wave  impact  on 
the  human  eye.  We  simulated  the  propagation  of  the  blast  wave  in  three-dimensional 
coordinates  to  obtain  the  transient  flow-fields  and  pressure  loading  on  the  human  eye. 
The  reflections  of  the  blast  wave  from  the  nasal  and  brow  ridges  amplify  the  pressure 
loading  and  produce  asymmetric  loading  on  the  eye.  We  predicted  possible  eye  injuries 
and  associated  injury  risks  using  the  calculated  intra-ocular  pressure  (IOP)  and  von-Mises 
stress  in  the  sclera.  The  sensitivity  study  of  the  model  to  material  properties  and 
anatomical  parameters  was  presented  and  discussed.  The  highest  injury  risk  from  blast 
loading  was  associated  with  cornea  abrasions  and  the  low  IOP  experienced  by  the  globe 
were  too  low  to  produce  an  appreciable  risk  for  globe  rupture.  The  maximum  IOP  was 
found  at  the  posterior  sclera,  which  may  increase  the  risk  for  delayed  injuries  to  the  optic 
nerve  head.  The  high  distortional  deformation  at  the  muscle  attachment  site  may  also 
cause  appreciable  risk  for  muscle  tearing  and  globe  dislocation. 
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